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\Q We investigate the dynamical properties of the 

transcriptional regulation of gene expression in 

^ the yeast Saccharomyces Cerevisiae within the 
^ framework of a synchronously and determin- 
1^ istically updated Boolean network model. By 
O means of a dynamically determinant subnet- 
^ work, we explore the robustness of transcrip- 
Q^tional regulation as a function of the type of 
' Boolean functions used in the model that mimic 
the influence of regulating agents on the tran- 
> scription level of a gene. We compare the re- 
^ suits obtained for the actual yeast network with 
^ those from two different model networks, one 
with similar in-degree distribution as the yeast 
CN and random otherwise, and another due to Bal- 
^ can et al. , where the global topology of the yeast 
network is reproduced faithfully. We, surpris- 
ingly, find that the first set of model networks 
^ better reproduce the results found with the ac- 



tual yeast network, even though the Balcan et 



d al. model networks are structurally more simi- 
lar to that of yeast. 



INTRODUCTION Recent advances in 
biotechnology allowed the accumulation 
of a vast amount of experimental data 
on intra-cellular processes, however, our 
knowledge on how a cell works re mains 
incomplete fLockh art & Winzelerl . 2000 : 
Barabasi & Olvai .20041] . The key compo- 
nent of the functional organization in a cell 
is the regulation of gene expression. By now, 
interacting gene pairs for several organisms 
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have been identified with 
age BBergmann et al 



significant cover- 
20031) . In particular 



the set of regulatory interactions identified in 
Saccharomyces Cerevisiae are believ ed to be 
close to complete llTeixeira et al. . I2OO6I1 . 



The activation/suppression dynamics in a cell 
evolves on a complex and inhomogeneous net- 
work of interactions. Therefore, the frame- 
work of graph theory serve as a power- 
ful mathematical tool for studying the reg- 
ulation of the gene expression on cellular 



level 


Albert & Barabasii 12002 


[ iBollobas', T998; 


Milo et al., 


2002; Colizza et al. 


, I2OO6; Newman, 


2001: 


Barabasi & Olvai. 


2004 


The topology of 



the graph describing the gene regulatory net- 
work (GRN) is far from being random and has 
been studied for sever al organisms, in partic- 



ular the budding yeast llGuelzim & et al 



Nicholas & et al.L 120041: iBergmann et al 



2OO2L 



20031. 



Deterministically and synchronously updated 
Boolean networks have been used widely as 
a mod el for r e gulato r y dynamics [ Kauffman, 



1969; lAldanaL I2003L iBalcan & ErzanL Iwoj 
I2OO6II . In this model, the expression levels 
of genes are discretized to take values or 
1 at each time ste p. Although i t is a major 

this 



2007D, 



oversimplification jNorrell et al 
approach has proven valuable in th e con- 
text of gene regula tion ^Mendoza et al.', 1999; 
Espinosa-Soto et all 12004: Albert & Othmer,. 



200311 . 



The network topology of the yeast's GRN is 
now believed to be unveiled to a large extent. 
However the nature of interactions, i.e., the 
rules that govern the dynamics, are not known 
in comparable detail. Accordingly, a statistical 
approach involving randomly assigned functions 
is relevant. Several classes of such functions 
have been investigated in the literature. The 



unbiased choice is to pick random Boolean func- 
tions. On the other hand, it has been claimed 
that experimental data is consistent with a sub- 
set of Boolean functions where one of the output 
is fixed for a particular value of one of the inp uts 



(canalizing functions) llHarris & et alj. 1200 



nputs 
m. It 



has also been suggested that a subset of canal- 
izing functions (nested canalizing functions) is 
more app ropriate for gene regu lation dynamics 
on yeast IIKauffman et alJ . 120031] . A more recent 
study finds that two subclasses of the nested 
canalizin g functions are actually dominant in 
the yeast IINikolaiewa et al.l . 120061] . 

The computational bottleneck in the analysis 
of Boolean network dynamics is the fact of that 
number of states increases exponentially with 
system size. This makes an exhaustive enumer- 
ation prohibitive, even if, in most cases, a frac- 
tion of the nodes can be left outside the anal- 
ysis due to their irrelevance to the dynamics 
by virtue of eith er the top olog y or the choic e of 
the function set I Socolar & KauffmanL 2003l] . In 
this paper, we determine and use a strongly con- 
nected subset of the genes that dictates the net- 
work's dynamical character and use a statistical 
approach to identify its robustness. 

The paper is organized as follows. In the 
Method section, we present the yeast's gene reg- 
ulation network, describe the employed Boolean 
dynamics, define the function classes used for 
setting the rules of the dynamics, propose a dy- 
namically relevant subnetwork and the model 
networks used for comparison with yeast. In the 
Result section, we present and analysis of the 
yeast's GRN dynamics, in particular exploring 
the robustness of the network under small per- 
turbations, comparing the results for different 
types of functions and for the actual vs. model 
network topologies. We discuss our findings in 
the last section. 



METHOD & MODELS Transcriptional reg- 
ulation of gene expression in a cell operates 
through transcription factors (TFs). These pro- 
teins bind the DNA on "promoter regions" (PRs) 
that act as the regulation centers of each gene. 
The details of this interaction can be very com- 
plex. In our study, as in past studies in the lit- 
erature, we assume that effect of the TFs that 
regulate a certain gene can be summarized in 
a Boolean function whose inputs represent the 
presence or the absence of TFs and the output 



determines whether the gene is activated or in- 
hibited for the given expression profile of the TF 
genes. 

The regulation dynamics evolves on a directed 
graph, whose nodes are the genes and a directed 
edge from A to B indicates that the product of 
A regulates B. The corresponding network for 
Saccharomyces Cerevisiae can be retrieve d from 
YEASTRACT database iTeixeira et al.L [2006] 
(www.yeastract.com). In order to be able to 
compare our results with past studies, we here 
consider an earlier version (2005) of the net- 
work including 4252 genes (with 146 TFs) with 
12541 interactions. As explained below, we 
also consider two model networks, one with a 
similar in-degree distribution as the yeast net- 
work above and random otherwise, and another 
with a topology highly similar to that of yeast, 
which emerg es from a n ull-model proposed ear- 
lier [ Balcan et al.Ll2007l] . 

The Boolean regulation dynamics on these 
networks is investigated by means of a syn- 
chronous and deterministic update of the net- 
work state as follows: Each node (gene) i has 
a state ai{t) at a particular time t where cri{t) 
is either 1 (on) or (off). The network state 
S{t) is the set of individual node states: S{t) = 
{cri(t), cj2(t), .., crAr(t)}. cri{t + 1) is determined by 
the Boolean function Bi assigned to i, which is 
a function of the states of the neighbor nodes 
connected to i by incoming edges. We used four 
types of random function classes found in the lit- 
erature as described below. 




Figure 1: The dynamically relevant 
(sub)network obtained after recursively pruning 
the (round) nodes with either zero out-degree or 
zero in-degree. 



1- Simple Random Function, RF: The 

rule table is constructed by setting the output 
for each input combination to 1 with probabilty 
p and otherwise, independent of the input. 

2- Canalizing Random Function, CF: 

A subclass of the RFs, that has at least 
one canalizing input variable whose canalizing 



value d etermin es the output IIKauffmanL Il993 
[kauffman et al.,. 20031 . 



a. 



(1) 



Bi{ai^i, ..,Sj, ..,ai^kin) / Sj 



where ]th in-neighbor is the canalizing node 
with Sj as the canalizing value and Sj as the 
canalizing output. Again, the output is de- 
termined through the parameter p. When 
the canalization condition is not satisfied, 
1, .., sJ, .., cjj^fc-^) in Exps. [T]is considered to 
be a RF. 



3- Nested Canalizing Random Func- 
tion, NCF: Nested Canalizing or Hierarchi- 
cally Canalizing functions are believed to bet- 
ter mo del gene regu l ation in biological sys- 
tems IKauffman et al.l. l2003ll . They form a sub- 
class of CFs, where one defines a canalizing or- 
der to the input nodes and the output is deter- 
mined by the first node in its canalizing value: 



Si,l 
Si,2 

Si,j 

Si,ki^ CTi,! 7^ Si A o-i,2 7^ S2 A ... A ai^k,„ = Ski„ 

. Sij^ <Tj,i 7^ Si A aj,2 / S2 A ... A ai^k,„ / s^^^ 



7^ si A o-j o = S2 



7^ Si A C7j,2 7^ S2 A ... A cr^ 



(2) 



We modify t he original definition in 

BKaufFman et al.L 12003 1. for the sake of an 
unbiased comparison with the other cases, 
by determining the outputs {sJ with the 
parameter p as before. 

4- Special Subclasses of Nested Canal- 
izing Random Function, SNCF: Follow- 
ing Nikolejewa et al., one can represent the 



NCFs above in a "minima l logical expres- 
sion" llNikolaiewa et al.l.l2006tl : 



^lfiOK2 0(-OK. 



-i0^f,kj-)) 



(3) 

where O represents either AND or OR logical 
function, i.e. Q G {A,V} and cr® stands for a 
possible negation of cr, i.e. cr® G ia, a|. Upon in- 



vestigation of Harris et al. data iHarris & et al 



[2b02]3 they found that gene regulatory rules are 
mainly governed by two subclasses of NCF: 



^f,l A K^2 A (... A K® A a^kj-)) 



(4) 



and 



-flA 



'i,2 



A ... A 



)...)) (5) 



with 66.39% and 29.41% probability of occur- 
rence, respectively. For these two functions, p is 
not a free parameter and depends on the topol- 
ogy 

Once the network topology and the functions 
are fixed, the Boolean dynamics is character- 
ized by a set of limit cycles which are the at- 
tractors of the dynamics reached from differ- 
ent initial conditions. Since these are the re- 
gions of the state space where the dynamics con- 
verges to, one expects them to be biologically 
relevant. For example, they have been associ- 
ated with different phenotypes of the pla nt Ara- 



bidopsis thaliana r Mendoza et al.L Il999 1 when 



the involved genes are those that take part 
in cell differentiation. We have investigated 
and compared the number, cycle-length, tran- 
sient length, and the basin of attraction of the 
attractors in each case. These results will 



be presented else where llTugrul & Kabakcioglu . 
expected in 2009tl . Here, we deal with another 
dynamical property, the robustness of the at- 
tractors to perturbation. We use the following 
arguments to measure th e robustness as given 
by Aldana IIAldanai 1200311 . Consider two copies 
of a network at states S{t) and S (t). Their Ham- 
ming Distance HD{t) is the number of nodes 
that differ between the two: 



N 



HD{t) = Y,\<^^{t)-<T[{t) 



(6) 



2=1 



Let x{t) = 1 — -^-S^, where N is the number of 



^private communication 



2 



• — • Simple Random Func. 




Figure 2: Robustness of yeast's GRN for all types of functions. For each p value, robustness was com- 
puted with averaging over 1000 random initial conditions of each 10 realization. Also, the Derrida's 
Exp., s = 2p{l - p){kin) was drawn. 



the nodes in the network. One defines the ro 
bustness s of the network as 

dx(t + 1) 



lim 

»1~ ,t— »oo 



dx{t) 



(7) 



The system is robust against perturbations (or- 
dered) if s < 1, whereas it is highly sen- 
sitive (chaotic) otherwise . It w as suggested 
by Kauffman jKauffman . I1993I1 that the ge- 
netic regulatory networks function at the edge 
of chaos, where s ~ 1. The quantity s can 
be estimated analytically for the RF case un- 
der an annealed approximation, both for ran- 
dom IDerrida & Pomeau, 1986] and power-law 
networks jAldanai 120031) . Derrida's result for 
random networks is 



2p{l-p){kin) 



(8) 



where p is, again, the unbiased probability that 
a binary function assigned to a node returns 1. 
Note that, by symmetry, s{p) = s{l—p). We mea- 
sure s numerically as a function of p € [0, 0.5], by 
examining the deviation of the two copies which 
are initially only slightly perturbed. For a net- 
work with nodes, the deviation is measured 
within a time window of 2N steps. 



As long as one is interested in the net- 
work characteristics such as attractor statistics 
or robustness, simulating the dynamics of the 
whole network is extremely inefficient. The 
reason is that, given the topology, some of 
the nodes make no contr ibution to such statis- 
tics I Socolar & KauffmanL [20031. For example, a 
node with zero in-degree remains at a fix state 
all times. Similarly, a node with zero out-degree 
simply follows the input and does not give any 
feedback. Same statements apply to those nodes 
which lose all incoming or outgoing edges after 
a round of pruning such nodes. Therefore, we 
focus on the dynamically relevant (sub)network 
(DRN) which is found by recursively pruning all 
the nodes with zero in-degree or zero out-degree 
(see Fig. [H). This subnetwork is typically much 
smaller than the original, allowing one to run 
time-efficient simulations. 

We find that the in-degree distribution of the 
DRN for the yeast regulatory network is expo- 
nential with an exponent a = 0.38, similar to the 
full yeast network. In addition, we generated 
two ensembles of 100 model networks for com- 
parison. The first ensemble is a set of randomly 




Figure 3: Robustness of in-EXP model networks with RF, CF and NCF type regulatory functions. The 
SNCF function results in s = 0.77 ± 0.05, whereas for the same p value yeast network has s = 0.78. 



connected networks of the same size = 82 and 
the same exponent in-degree distribution as the 
yeast DRN, so named in-EXP model. The second 
ensemble is generated by using a recently pro- 
posed null-model by Balcan et al. BBalcan et al. . 
200711 which successfully reproduces many topo- 
logical features of the yeast's GRN. An inter- 
esting observation is that, the second ensemble 
which preserves a number of topological signa- 
tures found in yeast, yielded dynamically rele- 
vant subnetworks with an average size of 36±15, 
i.e., significanly smaller than that of the yeast 
DRN. 



RESULTS On the GRN of the yeast, we cal- 
culated the robustness of the network dynamics 
for each function type discussed above as a func- 
tion of p. We chose p € {0.00, 0.01, .., 0.50}, ex- 
cept for SNCF case, where the value of p is fixed 
by the network topology to the value p = 0.27. 
For each case, we performed statistics over 10 
independent function assignments and ran the 
dynamics for 2N time steps starting from 1000 
different initial conditions (and their perturba- 
tions, in parallel). In all cases, the lengths of 
the simulations were sufficient for the system 
to reach an attractor. Fig.© shows the aver- 



age robustness found for each function type. We 
find that when a random RF function is associ- 
ated with each gene's transcription, the systems 
switches from an ordered phase to a chaotic 
phase around p = 0.22, consistent with Derrida's 
analytical result in Eq.®. CFs always result in 
an ordered system, except when p = 0.5, where 
the yeast's GRN appears to be at the edge of 
chaos. NCF and SNCFs which have been sug- 
gested to better represent gene regulation dy- 
namics are strictly ordered for all p values. 

For comparison, we repeated the same anal- 
ysis on in-EXP and Balcan et al. models. 100 
different networks were created from each set 
in order to reduce fluctuations due to structural 
deviations from sample to sample. The average 
robustness obtained for in-EXP and Balcan et 
al. model are compared with the corresponding 
data obtained from the yeast's GRN in Fig.® 
and Fig.®, respectively. We find that the in- 
EXP model networks show similar robustness 
profiles as the yeast's GRN in all cases, whereas 
Balcan et al. model, although it glob ally appears 



to capture the network structure [ Balcan et al 



,2007.1 . shows a significant deviation from the 
yeast in its dynamics. 
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Figure 4: Robustness of Balcan et al. model networks with RF, CF and NCF type regulatory func- 
tions. The SNCF function results in s = 0.83 ± 0.08, whereas for the same p value yeast network has 
s = 0.78. 



DISCUSSION We calculated the robustness 
of the yeast's transcriptional regulatory network 
within the framework of Boolean networks and 
as a function of the gene activation probability 
p. Under different assumptions on the func- 
tion class that governs the regulation process, 
we find that the network may show an order- 
chaos transition with changing p, may reach the 
edge of chaos at p = 0.5, or may stay robust for 
all p values. Our results point to the fact that, 
the activation probability by itself is not suffi- 
cient to determine the robustness of the Boolean 
networks; the functional category of the update 
rules also matter. As future experiments more 
precisely unveil activation/inhibition relations 
between genes in the yeast organism, proper 
choices for p and the function class shall become 
apparent. The strong dependence of the robust- 
ness to the function type and p may entail that 
both have been optimized throughout evolution- 
ary time scales to their present-time values. 
Our findings may then help addres s the "edge 
of chaos" hypothesis of Kauffman rtKauffmanl 



199311 



We furthermore compared our results on the 
yeast network with those obtained from two 



models which produce statistically similar net- 
work topologies. We found to our surprise that 
among the two models, Balcan et al. model 
which better reproduces a set of global topologi- 
cal features shows significantly larger deviation 
from the yeast's network in its robustness. This 
discrepancy should stem from certain structural 
features that are not captured by the global 
topological signatures considered in past stud- 
ies. One such difference we observe is the much 
smaller average dynamical core size of Balcan 
et al. model networks. This and other possible 
structural sources for the observed difference in 
dynamics should be further investigated. 
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